home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / DCLAP 6d / dclap6d / corelib / ncbimath.c < prev    next >
Text File  |  1996-07-05  |  18KB  |  754 lines
  1. /*   ncbimath.c
  2. * ===========================================================================
  3. *
  4. *                            PUBLIC DOMAIN NOTICE
  5. *               National Center for Biotechnology Information
  6. *
  7. *  This software/database is a "United States Government Work" under the
  8. *  terms of the United States Copyright Act.  It was written as part of
  9. *  the author's official duties as a United States Government employee and
  10. *  thus cannot be copyrighted.  This software/database is freely available
  11. *  to the public for use. The National Library of Medicine and the U.S.
  12. *  Government have not placed any restriction on its use or reproduction.
  13. *
  14. *  Although all reasonable efforts have been taken to ensure the accuracy
  15. *  and reliability of the software and data, the NLM and the U.S.
  16. *  Government do not and cannot warrant the performance or results that
  17. *  may be obtained by using this software or data. The NLM and the U.S.
  18. *  Government disclaim all warranties, express or implied, including
  19. *  warranties of performance, merchantability or fitness for any particular
  20. *  purpose.
  21. *
  22. *  Please cite the author in any work or product based on this material.
  23. *
  24. * ===========================================================================
  25. *
  26. * File Name:  ncbimath.c
  27. *
  28. * Author:  Gish, Kans, Ostell, Schuler
  29. *
  30. * Version Creation Date:   10/23/91
  31. *
  32. * $Revision: 2.7 $
  33. *
  34. * File Description:
  35. *       portable math functions
  36. *
  37. * Modifications:
  38. * --------------------------------------------------------------------------
  39. * Date     Name        Description of modification
  40. * -------  ----------  -----------------------------------------------------
  41. * 04-15-93 Schuler     Changed _cdecl to LIBCALL
  42. * 12-22-93 Schuler     Converted ERRPOST((...)) to ErrPostEx(...)
  43. *
  44. * ==========================================================================
  45. */
  46.  
  47. #undef  THIS_MODULE
  48. #define THIS_MODULE g_corelib
  49. #define THI_FILE _this_file
  50.  
  51. #include <ncbi.h>
  52. #include <ncbiwin.h>
  53.  
  54. extern char * g_corelib;
  55. static char * _this_file = __FILE__;
  56.  
  57. /*
  58.     Nlm_Expm1(x)
  59.     Return values accurate to approx. 16 digits for the quantity exp(x)-1
  60.     for all x.
  61. */
  62. double LIBCALL
  63. Nlm_Expm1(x)
  64.     register double    x;
  65. {
  66.     register double    absx;
  67.  
  68.     if ((absx = ABS(x)) > .33)
  69.         return exp(x) - 1.;
  70.  
  71.     if (absx < 1.e-16)
  72.         return x;
  73.  
  74.     return x * (1. + x *
  75.                     (0.5 + x * (1./6. + x *
  76.                         (1./24. + x * (1./120. + x *
  77.                             (1./720. + x * (1./5040. + x *
  78.                                 (1./40320. + x * (1./362880. + x *
  79.                                     (1./3628800. + x * (1./39916800. + x *
  80.                                         (1./479001600. + x/6227020800.)
  81.                                     ))
  82.                                 ))
  83.                             ))
  84.                         ))
  85.                     )));
  86. }
  87.  
  88. /*
  89.     Nlm_Log1p(x)
  90.     Return accurate values for the quantity log(x+1) for all x > -1.
  91. */
  92. double LIBCALL
  93. Nlm_Log1p(x)
  94.     register double    x;
  95. {
  96.     register int    i;
  97.     register double    sum, y;
  98.  
  99.     if (ABS(x) >= 0.2)
  100.         return log(x+1.);
  101.  
  102.     for (i=0, sum=0., y=x; ; ) {
  103.         sum += y/++i;
  104.         if (y < DBL_EPSILON)
  105.             break;
  106.         y *= x;
  107.         sum -= y/++i;
  108.         if (y < DBL_EPSILON)
  109.             break;
  110.         y *= x;
  111.     }
  112.     return sum;
  113. }
  114.  
  115.  
  116. /*
  117. ** Special thanks to Dr. John ``Gammas Galore'' Spouge for deriving the
  118. ** method for calculating the gamma coefficients and their use.
  119. ** (See the #ifdef-ed program included at the end of this file).
  120. **/
  121.  
  122. /*
  123.     For discussion of the Gamma function, see "Numerical Recipes in C",
  124.     Press et al. (1988) pages 167-169.
  125. */
  126.  
  127. static double NEAR general_lngamma PROTO((double x,int n));
  128.  
  129. static double _default_gamma_coef [] = {
  130.      4.694580336184385e+04,
  131.     -1.560605207784446e+05,
  132.      2.065049568014106e+05,
  133.     -1.388934775095388e+05,
  134.      5.031796415085709e+04,
  135.     -9.601592329182778e+03,
  136.      8.785855930895250e+02,
  137.     -3.155153906098611e+01,
  138.      2.908143421162229e-01,
  139.     -2.319827630494973e-04,
  140.      1.251639670050933e-10
  141.      };
  142. static double    PNTR gamma_coef = _default_gamma_coef;
  143. static unsigned    gamma_dim = DIM(_default_gamma_coef);
  144. static double    xgamma_dim = DIM(_default_gamma_coef);
  145.  
  146. void LIBCALL
  147. Nlm_GammaCoeffSet(double PNTR cof, unsigned dim) /* changes gamma coeffs */
  148. {
  149.     if (dim < 3 || dim > 100) /* sanity check */
  150.         return;
  151.     gamma_coef = cof;
  152.     xgamma_dim = gamma_dim = dim;
  153. }
  154.  
  155.  
  156. static double NEAR
  157. general_lngamma(x, order)    /* nth derivative of ln[gamma(x)] */
  158.     double    x;        /* 10-digit accuracy achieved only for 1 <= x */
  159.     int        order;    /* order of the derivative, 0..POLYGAMMA_ORDER_MAX */
  160. {
  161.     int        i;
  162.     double    xx, tx;
  163.     double    y[POLYGAMMA_ORDER_MAX+1];
  164.     register double    tmp, value, PNTR coef;
  165.  
  166.     xx = x - 1.;  /* normalize from gamma(x + 1) to xx! */
  167.     tx = xx + xgamma_dim;
  168.     for (i = 0; i <= order; ++i) {
  169.         tmp = tx;
  170.         /* sum the least significant terms first */
  171.         coef = &gamma_coef[gamma_dim];
  172.         if (i == 0) {
  173.             value = *--coef / tmp;
  174.             while (coef > gamma_coef)
  175.                 value += *--coef / --tmp;
  176.         }
  177.         else {
  178.             value = *--coef / Nlm_Powi(tmp, i + 1);
  179.             while (coef > gamma_coef)
  180.                 value += *--coef / Nlm_Powi(--tmp, i + 1);
  181.             tmp = Nlm_Factorial(i);
  182.             value *= (i%2 == 0 ? tmp : -tmp);
  183.         }
  184.         y[i] = value;
  185.     }
  186.     ++y[0];
  187. #if !defined(DCLAP) || !defined(COMP_CWI)
  188.     value = Nlm_LogDerivative(order, y);
  189. #endif
  190.     tmp = tx + 0.5;
  191.     switch (order) {
  192.     case 0:
  193.         value += ((NCBIMATH_LNPI+NCBIMATH_LN2) / 2.)
  194.                     + (xx + 0.5) * log(tmp) - tmp;
  195.         break;
  196.     case 1:
  197.         value += log(tmp) - xgamma_dim / tmp;
  198.         break;
  199.     case 2:
  200.         value += (tmp + xgamma_dim) / (tmp * tmp);
  201.         break;
  202.     case 3:
  203.         value -= (1. + 2.*xgamma_dim / tmp) / (tmp * tmp);
  204.         break;
  205.     case 4:
  206.         value += 2. * (1. + 3.*xgamma_dim / tmp) / (tmp * tmp * tmp);
  207.         break;
  208.     default:
  209.         tmp = Nlm_Factorial(order - 2) * Nlm_Powi(tmp, 1 - order)
  210.                 * (1. + (order - 1) * xgamma_dim / tmp);
  211.         if (order % 2 == 0)
  212.             value += tmp;
  213.         else
  214.             value -= tmp;
  215.         break;
  216.     }
  217.     return value;
  218. }
  219.  
  220.  
  221. double LIBCALL
  222. Nlm_PolyGamma(x, order) /* ln(ABS[gamma(x)]) - 10 digits of accuracy */
  223.     double    x;      /* and derivatives */
  224.     int        order;    /* order of the derivative */
  225. /* order = 0, 1, 2, ...  ln(gamma), digamma, trigamma, ... */
  226. /* CAUTION:  the value of order is one less than the suggested "di" and
  227. "tri" prefixes of digamma, trigamma, etc.  In other words, the value
  228. of order is truly the order of the derivative.  */
  229. {
  230.     int        k;
  231.     register double    value, tmp;
  232.     double    y[POLYGAMMA_ORDER_MAX+1], sx;
  233.  
  234.     if (order < 0 || order > POLYGAMMA_ORDER_MAX) {
  235.         ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"PolyGamma: unsupported derivative order");
  236.         /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "unsupported derivative order"));**/
  237.         return HUGE_VAL;
  238.     }
  239.  
  240.     if (x >= 1.)
  241.         return general_lngamma(x, order);
  242.  
  243.     if (x < 0.) {
  244.         value = general_lngamma(1. - x, order);
  245.         value = ((order - 1) % 2 == 0 ? value : -value);
  246.         if (order == 0) {
  247.             sx = sin(NCBIMATH_PI * x);
  248.             sx = ABS(sx);
  249.             if ( (x < -0.1 && (ceil(x) == x || sx < 2.*DBL_EPSILON))
  250.                     || sx == 0.) {
  251.                 ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"PolyGamma: log(0)");
  252.                 /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "log(0)"));**/
  253.                 return HUGE_VAL;
  254.             }
  255.             value += NCBIMATH_LNPI - log(sx);
  256.         }
  257.         else {
  258.             y[0] = sin(x *= NCBIMATH_PI);
  259.             tmp = 1.;
  260.             for (k = 1; k <= order; k++) {
  261.                 tmp *= NCBIMATH_PI;
  262.                 y[k] = tmp * sin(x += (NCBIMATH_PI/2.));
  263.             }
  264. #if !defined(DCLAP) || !defined(COMP_CWI)
  265.             value -= Nlm_LogDerivative( order, y);
  266. #endif
  267.         }
  268.     }
  269.     else {
  270.         value = general_lngamma(1. + x, order);
  271.         if (order == 0) {
  272.             if (x == 0.) {
  273.                 ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"PolyGamma: log(0)");
  274.                 /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "log(0)"));**/
  275.                 return HUGE_VAL;
  276.             }
  277.             value -= log(x);
  278.         }
  279.         else {
  280.             tmp = Nlm_Factorial(order - 1) * Nlm_Powi(x,  -order);
  281.             value += (order % 2 == 0 ? tmp : - tmp);
  282.         }
  283.     }
  284.     return value;
  285. }
  286.  
  287. double LIBCALL Nlm_LogDerivative(order, u) /* nth derivative of ln(u) */
  288.     int        order;        /* order of the derivative */
  289.     double    PNTR u;        /* values of u, u', u", etc. */
  290. {
  291.     int        i;
  292.     double    y[LOGDERIV_ORDER_MAX+1];
  293.     register double    value, tmp;
  294.  
  295.     if (order < 0 || order > LOGDERIV_ORDER_MAX) {
  296. InvalidOrder:
  297.         ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"LogDerivative: unsupported derivative order");
  298.         /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "unsupported derivative order"));**/
  299.         return HUGE_VAL;
  300.     }
  301.  
  302.     if (order > 0 && u[0] == 0.) {
  303.         ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"LogDerivative: divide by 0");
  304.         /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "divide by 0"));**/
  305.         return HUGE_VAL;
  306.     }
  307.     for (i = 1; i <= order; i++)
  308.         y[i] = u[i] / u[0];
  309.  
  310.     switch (order) {
  311.     case 0:
  312.         if (u[0] > 0.)
  313.             value = log(u[0]);
  314.         else {
  315.             ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"LogDerivative: log(x<=0)");
  316.             /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "log(x<=0)"));**/
  317.             return HUGE_VAL;
  318.         }
  319.         break;
  320.     case 1:
  321.         value = y[1];
  322.         break;
  323.     case 2:
  324.         value = y[2] - y[1] * y[1];
  325.         break;
  326.     case 3:
  327.         value = y[3] - 3. * y[2] * y[1] + 2. * y[1] * y[1] * y[1];
  328.         break;
  329.     case 4:
  330.         value = y[4] - 4. * y[3] * y[1] - 3. * y[2] * y[2]
  331.                           + 12. * y[2] * (tmp = y[1] * y[1]);
  332.         value -= 6. * tmp * tmp;
  333.         break;
  334.     default:
  335.         goto InvalidOrder;
  336.     }
  337.     return value;
  338. }
  339.  
  340.  
  341.  
  342. double LIBCALL
  343. Nlm_Gamma(x)        /* ABS[gamma(x)] - 10 digits of accuracy */
  344.     double     x;
  345. {
  346.     return exp(Nlm_PolyGamma(x, 0));
  347. }
  348.  
  349. double LIBCALL
  350. Nlm_LnGamma(x)        /* ln(ABS[gamma(x)]) - 10 digits of accuracy */
  351.     double     x;
  352. {
  353.     return Nlm_PolyGamma(x, 0);
  354. }
  355.  
  356. double LIBCALL
  357. Nlm_DiGamma(x)    /* digamma, 1st order derivative of log(gamma) */
  358.     double    x;
  359. {
  360.     return Nlm_PolyGamma(x, 1);
  361. }
  362.  
  363. double LIBCALL
  364. Nlm_TriGamma(x)    /* trigamma, 2nd order derivative of log(gamma) */
  365.     double    x;
  366. {
  367.     return Nlm_PolyGamma(x, 2);
  368. }
  369.  
  370.  
  371. #ifdef foo
  372. /* A program to calculate the gamma coefficients based on a method
  373.    by John Spouge.
  374.    Cut this program out, save it in a separate file, and compile.
  375.    Be sure to link with a math library.
  376. */
  377. /*
  378.     a[n] = ((gamma+0.5-n)^(n-0.5)) * exp(gamma+0.5-n) *
  379.         ((-1)^(n-1) / (n-1)!) * (1/sqrt(2*Pi))
  380. */
  381. #include <stdio.h>
  382. #include <math.h>
  383.  
  384. main(ac, av)
  385.     int ac;
  386.     char **av;
  387.  
  388. {
  389.     int    i, j, cnt;
  390.     double    a, x, y, z, ifact = 1.;
  391.  
  392.     if (ac != 2 || sscanf(av[1], "%d", &cnt) != 1)
  393.         exit(1);
  394.  
  395.     for (i=0; i<cnt; ++i) {
  396.         x = cnt - (i + 0.5);
  397.         y = log(x) * (i + 0.5) + x;
  398.         y = exp(y);
  399.         if (i > 1)
  400.             ifact *= i;
  401.         y /= ifact;
  402.         if (i%2 == 1)
  403.             y = -y;
  404.         printf("\t\t\t%.17lg", y);
  405.         if (i < cnt-1)
  406.             putchar(',');
  407.         putchar('\n');
  408.     }
  409. }
  410. #endif /* foo */
  411.  
  412.  
  413. #define FACTORIAL_PRECOMPUTED    36
  414.  
  415. double LIBCALL
  416. Nlm_Factorial(int n)
  417.  
  418. {
  419.     static double    precomputed[FACTORIAL_PRECOMPUTED]
  420.                 = { 1., 1., 2., 6., 24., 120., 720.};
  421.     static int    nlim = 6;
  422.  
  423.     if (n < 0)
  424.         return 0.0; /* Undefined! */
  425.  
  426.     if (n >= DIM(precomputed))
  427.         return exp(Nlm_LnGamma((double)(n+1)));
  428.  
  429.     while (nlim < n) {
  430.         ++nlim;
  431.         precomputed[nlim] = precomputed[nlim-1] * nlim;
  432.     }
  433.     return precomputed[n];
  434. }
  435.  
  436. /* Nlm_LnGammaInt(n) -- return log(Gamma(n)) for integral n */
  437. double LIBCALL
  438. Nlm_LnGammaInt(int n)
  439. {
  440.     static double    precomputed[FACTORIAL_PRECOMPUTED];
  441.     static int    nlim = 1; /* first two entries are 0 */
  442.  
  443.     if (n >= 0 && n < DIM(precomputed)) {
  444.         while (nlim < n) {
  445.             ++nlim;
  446.             precomputed[nlim] = log(Nlm_Factorial(nlim-1));
  447.         }
  448.         return precomputed[n];
  449.     }
  450.     return Nlm_LnGamma((double)n);
  451. }
  452.  
  453.  
  454.  
  455. /*
  456.     Combined Newton-Raphson and Bisection root-finder
  457.  
  458.     Original Function Name:  Inv_Xnrbis()
  459.     Author:  Dr. John Spouge
  460.     Location:  NCBI
  461.     Received:  July 16, 1991
  462. */
  463. #define F(x)  ((*f)(x)-y)
  464. #define DF(x) ((*df)(x))
  465. #define NRBIS_ITMAX 100
  466.  
  467. double LIBCALL
  468. Nlm_NRBis(double y, double (LIBCALL *f )PROTO ((double )), double (LIBCALL *df )PROTO ((double )), double p, double x, double q, double tol)        /* tolerance */
  469.  
  470. {
  471.     double    temp;    /* for swapping end-points if necessary */
  472.     double    dx;        /* present interval length */
  473.     double    dxold;    /* old interval length */
  474.     double    fx;        /* f(x)-y */
  475.     double    dfx;    /* Df(x) */
  476.     int    j;                       /* loop index */
  477.     double    fp, fq;
  478.  
  479. /* Preliminary checks for improper bracketing and end-point root. */
  480.     if ((fp = F(p)) == 0.)
  481.         return p;
  482.     if ((fq = F(q)) == 0.)
  483.         return q;
  484.     if ((fp > 0. && fq > 0.) || (fp < 0. && fq < 0.)) {
  485.         ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_INVAL,"NRBis: root not bracketed");
  486.         /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_INVAL, "root not bracketed"));**/
  487.         return HUGE_VAL;
  488.     }
  489. /* Swaps end-points if necessary to make F(p)<0.<F(q) */
  490.     if (fp > 0.) {
  491.         temp = p;
  492.         p = q;
  493.         q = temp;
  494.     }
  495. /* Set up the Bisection & Newton-Raphson iteration. */
  496.     if ((x-p) * (x-q) > 0.)
  497.         x = 0.5*(p+q);
  498.     dxold = dx = p-q;
  499.     for (j = 1; j <= NRBIS_ITMAX; ++j) {
  500.         fx = F(x);
  501.         if (fx == 0.) /* Check for termination. */
  502.             return x;
  503.         if (fx < 0.)
  504.             p = x;
  505.         else
  506.             q = x;
  507.         dfx = DF(x);
  508. /* Check: root out of bounds or bisection faster than Newton-Raphson? */
  509.         if ((dfx*(x-p)-fx)*(dfx*(x-q)-fx) >= 0. || 2.*ABS(fx) > ABS(dfx*dx)) {
  510.             dx = dxold;               /* Bisect */
  511.             dxold = 0.5*(p-q);
  512.             x = 0.5*(p+q);
  513.             if (ABS(dxold) <= tol)
  514.                 return x;
  515.         } else {
  516.             dx = dxold;               /* Newton-Raphson */
  517.             dxold = fx/dfx;
  518.             x -= dxold;
  519.             if (ABS(dxold) < tol && F(x-SIGN(dxold)*tol)*fx < 0.)
  520.                 return x;
  521.         }
  522.     }
  523.     
  524.     ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_ITER,"NRBis: iterations > NRBIS_ITMAX");
  525.     /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_ITER, "iterations > NRBIS_ITMAX"));**/
  526.     return HUGE_VAL;
  527. }
  528. #undef F /* clean-up */
  529. #undef DF /* clean-up */
  530.  
  531.  
  532.  
  533. /*
  534.     Romberg numerical integrator
  535.  
  536.     Author:  Dr. John Spouge, NCBI
  537.     Received:  July 17, 1992
  538.     Reference:
  539.         Francis Scheid (1968)
  540.         Schaum's Outline Series
  541.         Numerical Analysis, p. 115
  542.         McGraw-Hill Book Company, New York
  543. */
  544. #define F(x)  ((*f)((x), fargs))
  545. #define ROMBERG_ITMAX 20
  546.  
  547. double LIBCALL
  548. Nlm_RombergIntegrate(double (LIBCALL *f) (double,Nlm_VoidPtr), Nlm_VoidPtr fargs, double p, double q, double eps, int epsit, int itmin)
  549.  
  550. {
  551.     double    romb[ROMBERG_ITMAX];    /* present list of Romberg values */
  552.     double    h;    /* mesh-size */
  553.     int        i, j, k, npts;
  554.     long    n;    /* 4^(error order in romb[i]) */
  555.     int        epsit_cnt = 0, epsck;
  556.     register double    y;
  557.     register double    x;
  558.     register double    sum;
  559.  
  560.     /* itmin = min. no. of iterations to perform */
  561.     itmin = MAX(1, itmin);
  562.     itmin = MIN(itmin, ROMBERG_ITMAX-1);
  563.  
  564.     /* epsit = min. no. of consecutive iterations that must satisfy epsilon */
  565.     epsit = MAX(epsit, 1); /* default = 1 */
  566.     epsit = MIN(epsit, 3); /* if > 3, the problem needs more prior analysis */
  567.  
  568.     epsck = itmin - epsit;
  569.  
  570.     npts = 1;
  571.     h = q - p;
  572.     x = F(p);
  573.     if (ABS(x) == HUGE_VAL)
  574.         return x;
  575.     y = F(q);
  576.     if (ABS(y) == HUGE_VAL)
  577.         return y;
  578.     romb[0] = 0.5 * h * (x + y);    /* trapezoidal rule */
  579.     for (i = 1; i < ROMBERG_ITMAX; ++i, npts *= 2, h *= 0.5) {
  580.         sum = 0.;    /* sum of ordinates for */
  581.         /* x = p+0.5*h, p+1.5*h, ..., q-0.5*h */
  582.         for (k = 0, x = p+0.5*h; k < npts; k++, x += h) {
  583.             y = F(x);
  584.             if (ABS(y) == HUGE_VAL)
  585.                 return y;
  586.             sum += y;
  587.         }
  588.         romb[i] = 0.5 * (romb[i-1] + h*sum); /* new trapezoidal estimate */
  589.  
  590.         /* Update Romberg array with new column */
  591.         for (n = 4, j = i-1; j >= 0; n *= 4, --j)
  592.             romb[j] = (n*romb[j+1] - romb[j]) / (n-1);
  593.  
  594.         if (i > epsck) {
  595.             if (ABS(romb[1] - romb[0]) > eps * ABS(romb[0])) {
  596.                 epsit_cnt = 0;
  597.                 continue;
  598.             }
  599.             ++epsit_cnt;
  600.             if (i >= itmin && epsit_cnt >= epsit)
  601.                 return romb[0];
  602.         }
  603.     }
  604.     ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_ITER,"RombergIntegrate: iterations > ROMBERG_ITMAX");
  605.     /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_ITER, "iterations > ROMBERG_ITMAX"));**/
  606.     return HUGE_VAL;
  607. }
  608.  
  609. /*
  610.     Nlm_Gcd(a, b)
  611.  
  612.     Return the greatest common divisor of a and b.
  613.  
  614.     Adapted 8-15-90 by WRG from code by S. Altschul.
  615. */
  616. long LIBCALL
  617. Nlm_Gcd(register long a, register long b)
  618. {
  619.     register long    c;
  620.  
  621.     b = ABS(b);
  622.     if (b > a)
  623.         c=a, a=b, b=c;
  624.  
  625.     while (b != 0) {
  626.         c = a%b;
  627.         a = b;
  628.         b = c;
  629.     }
  630.     return a;
  631. }
  632.  
  633. /* Round a floating point number to the nearest integer */
  634. long LIBCALL
  635. Nlm_Nint(register double x)    /* argument */
  636. {
  637.     x += (x >= 0. ? 0.5 : -0.5);
  638.     return (long)x;
  639. }
  640.  
  641. /*
  642. integer power function
  643.  
  644. Original submission by John Spouge, 6/25/90
  645. Added to shared library by WRG
  646. */
  647. double LIBCALL
  648. Nlm_Powi(register double x, register int n)    /* power */
  649. {
  650.     register double    y;
  651.  
  652.     if (n == 0)
  653.         return 1.;
  654.  
  655.     if (x == 0.) {
  656.         if (n < 0) {
  657.             ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"Powi: divide by 0");
  658.             /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "divide by 0"));**/
  659.             return HUGE_VAL;
  660.         }
  661.         return 0.;
  662.     }
  663.  
  664.     if (n < 0) {
  665.         x = 1./x;
  666.         n = -n;
  667.     }
  668.  
  669.     while (n > 1) {
  670.         if (n&1) {
  671.             y = x;
  672.             goto Loop2;
  673.         }
  674.         n /= 2;
  675.         x *= x;
  676.     }
  677.     return x;
  678.  
  679. Loop2:
  680.     n /= 2;
  681.     x *= x;
  682.     while (n > 1) {
  683.         if (n&1)
  684.             y *= x;
  685.         n /= 2;
  686.         x *= x;
  687.     }
  688.     return y * x;
  689. }
  690.  
  691. /*
  692.     Additive random number generator
  693.  
  694.     Modelled after "Algorithm A" in
  695.     Knuth, D. E. (1981). The art of computer programming, volume 2, page 27.
  696.  
  697.     7/26/90 WRG
  698. */
  699.  
  700. static long    state[33] = {
  701.     (long)0xd53f1852,  (long)0xdfc78b83,  (long)0x4f256096,  (long)0xe643df7,
  702.     (long)0x82c359bf,  (long)0xc7794dfa,  (long)0xd5e9ffaa,  (long)0x2c8cb64a,
  703.     (long)0x2f07b334,  (long)0xad5a7eb5,  (long)0x96dc0cde,  (long)0x6fc24589,
  704.     (long)0xa5853646,  (long)0xe71576e2,  (long)0xdae30df,  (long)0xb09ce711,
  705.     (long)0x5e56ef87,  (long)0x4b4b0082,  (long)0x6f4f340e,  (long)0xc5bb17e8,
  706.     (long)0xd788d765,  (long)0x67498087,  (long)0x9d7aba26,  (long)0x261351d4,
  707.     (long)0x411ee7ea,  (long)0x393a263,  (long)0x2c5a5835,  (long)0xc115fcd8,
  708.     (long)0x25e9132c,  (long)0xd0c6e906,  (long)0xc2bc5b2d,  (long)0x6c065c98,
  709.     (long)0x6e37bd55 };
  710.  
  711. #define r_off    12
  712.  
  713. static long    *rJ = &state[r_off],
  714.             *rK = &state[DIM(state)-1];
  715.  
  716. void LIBCALL
  717. Nlm_RandomSeed(long x)
  718. {
  719.     register int    i;
  720.  
  721.     state[0] = x;
  722.     /* linear congruential initializer */
  723.     for (i=1; i<DIM(state); ++i) {
  724.         state[i] = 1103515245*state[i-1] + 12345;
  725.     }
  726.  
  727.     rJ = &state[r_off];
  728.     rK = &state[DIM(state)-1];
  729.  
  730.     for (i=0; i<10*DIM(state); ++i)
  731.         (void) Nlm_RandomNum();
  732. }
  733.  
  734.  
  735. /*
  736.     Nlm_RandomNum --  return value in the range 0 <= x <= 2**31 - 1
  737. */
  738. long LIBCALL
  739. Nlm_RandomNum(void)
  740. {
  741.     register long    r;
  742.  
  743.     r = *rK;
  744.     r += *rJ--;
  745.     *rK-- = r;
  746.  
  747.     if (rK < state)
  748.         rK = &state[DIM(state)-1];
  749.     else
  750.         if (rJ < state)
  751.             rJ = &state[DIM(state)-1];
  752.     return (r>>1)&0x7fffffff; /* discard the least-random bit */
  753. }
  754.